1. Field of the Invention
This invention relates to a diffusion status prediction method for a diffused substance. The invention is designed to predict how a substance (for example, a radioactive substance or smoke), emitted from a diffusion source (for example, facilities using a radioactive substance, and a chimney) into the atmosphere, will diffuse in the atmosphere, thereby estimating the concentrations of the substance varying with time at different locations.
2. Description of the Related Art
An atmospheric diffusion status prediction method is under development which, if a radioactive substance is emitted to the outside by accident of facilities dealing with the radioactive substance, estimates the diffusion range of the radioactive substance and the concentrations of the radioactive substance at different locations, thereby predicting areas which may be endangered by the radioactive substance.
This diffusion status prediction method can be applied not only in predicting the diffusion status of the radioactive substance, but also in calculating the concentrations of a gaseous body (smoke) at different locations if the gaseous body emitted from a chimney of a factory diffuses into the atmosphere, as well as in analyzing the diffusion status of an air pollutant for environmental impact assessment.
To predict the diffusion status of the substance, emitted into the atmosphere, by computation, it is necessary to carry out the following two computations:
(1) Meteorology predictive computation
(2) Diffusion status predictive computation
The meteorology predictive computation (1) is carried out as follows: Partial differential equations, which analyze atmospheric phenomena, are solved based on meteorological observation data such as those from meteorological GPV (Grid Point Value) database in Japan and meteorological databases of foreign countries. The solutions give wind directions and wind speeds at many evaluation locations (grid point positions) at constant time intervals during the period of time from the initial stage of an event (e.g., emission of a radioactive substance to the outside) until a point in time which is a predetermined time after the occurrence of the event. That is, the meteorology predictive computation (1) finds an atmospheric condition representing wind field data at constant time intervals.
The diffusion status predictive computation (2) is carried out as follows: The concentration and properties of the diffused substance emitted, and the above wind field data are substituted into a diffusion equation for computing the diffused state of the substance (particles), thereby obtaining the concentrations of the diffused substance at the respective grid point positions at the respective time intervals.
<Explanation for the Meteorology Predictive Computation>
The outline of the meteorology predictive computation will be offered first. The meteorological observation data, for example, meteorological GPV data, are delivered online from the Meteorological Business Support Center every 12 hours. The meteorological GPV data show meteorological data (wind velocity vectors (wind direction, wind speed), atmospheric pressure, temperature, moisture content) at locations (called “parent grid point positions”) where a plurality of longitudinal virtual lines, which extend along a north-south direction on the earth's surface and which are spaced from each other in an east-west direction on the earth's surface by a prescribed distance (20 Km), cross a plurality of latitudinal virtual lines, which extend along the east-west direction on the earth's surface and which are spaced from each other in the north-south direction on the earth's surface by a prescribed distance (20 Km). In addition, the meteorological GPV data are delivered online as the meteorological data at the respective parent grid point positions which comprise a total of 51 hours of data obtained at 3-hour intervals, such as at the start of delivery, and 3 hours after, 6 hours after, and 9 hours after the start of delivery.
The above-mentioned meteorological data at the parent grid point positions, as the meteorological GPV data, involve a long distance of 20 Km between the parent grid point positions in terms of space, and a long interval of 3 hours in terms of time. Thus, the diffusion concentrations of the diffused substance cannot be computed only based on the gas status (wind direction, wind speed) data shown by the meteorological data at these parent grid point positions, namely, wind field data.
Hence, there is need to determine the atmospheric condition (wind direction, wind speed, etc.), which is detailed in terms of time as well as space, from the meteorological observation data rough spatially and rough timewise, by performing computations according to partial differential equations for analyzing atmospheric phenomena.
Thus, child grid point positions are set between parent grid point positions set in calculation areas to be calculated (specific areas preset in the earth's surface). The child grid point positions are arranged at the locations where a plurality of longitudinal virtual lines, which extend along the north-south direction on the earth's surface and which are spaced from each other in the east-west direction on the earth's surface by a constant distance (50 m), cross a plurality of latitudinal virtual lines, which extend along the east-wet direction on the earth's surface and which are spaced from each other in the north-south direction on the earth's surface by a constant distance (50 m).
The meteorological data at the child grid point positions and the parent grid point positions at predetermined time intervals (for example, 20-second intervals) from the start of computation are obtained by difference analysis computation according to partial differential equations, which analyze atmospheric phenomena. As the partial differential equations for analyzing atmospheric phenomena, there can be used fundamental equations for wind field analysis, which are represented by the RAMS (Regional Atmospheric Modeling System) code developed by Colorado State University and ATMET (a company in Colorado, USA).
The fundamental equations for wind field analysis, represented in the RAMS code, comprise equations of motion, an equation of thermal energy, an equation of moisture diffusion, and an equation of continuity, and are expressed as the following equations (1) to (6):                               Equations          ⁢                                           ⁢          of          ⁢                                           ⁢          motion                ⁢                                                                                                              ∂            u                                ∂            t                          =                                            -              u                        ⁢                                          ∂                u                                            ∂                x                                              -                      v            ⁢                                          ∂                u                                            ∂                y                                              -                      w            ⁢                                          ∂                u                                            ∂                z                                              -                      θ            ⁢                                          ∂                                  π                  ′                                                            ∂                x                                              +                      f            ⁢                                                   ⁢            v                    +                                                    ∂                                                                                               ∂                x                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    u                                                        ∂                    x                                                              )                                +                                    ∂                              ∂                y                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    u                                                        ∂                    y                                                              )                                +                                                    ∂                                                                                               ∂                z                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    u                                                        ∂                    z                                                              )                                                          (        1        )                                                      ∂            v                                ∂            t                          =                                            -              u                        ⁢                                          ∂                v                                            ∂                x                                              -                      v            ⁢                                          ∂                v                                            ∂                y                                              -                      w            ⁢                                          ∂                v                                            ∂                z                                              -                      θ            ⁢                                          ∂                                  π                  ′                                                            ∂                y                                              -                      f            ⁢                                                   ⁢            u                    +                                                    ∂                                                                                               ∂                x                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    v                                                        ∂                    x                                                              )                                +                                    ∂                              ∂                y                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    v                                                        ∂                    y                                                              )                                +                                                    ∂                                                                                               ∂                z                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    v                                                        ∂                    z                                                              )                                                          (        2        )                                                      ∂            w                                ∂            t                          =                                            -              u                        ⁢                                          ∂                w                                            ∂                x                                              -                      v            ⁢                                          ∂                w                                            ∂                y                                              -                      w            ⁢                                          ∂                w                                            ∂                z                                              -                      θ            ⁢                                          ∂                                  π                  ′                                                            ∂                z                                              -                                    g              ⁢                                                           ⁢                              θ                v                ′                                                    θ              0                                +                                                    ∂                                                                                               ∂                x                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    w                                                        ∂                    x                                                              )                                +                                                    ∂                                                                                               ∂                y                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    w                                                        ∂                    y                                                              )                                +                                                    ∂                                                                                               ∂                z                                      ⁢                          (                                                K                  m                                ⁢                                                      ∂                    w                                                        ∂                    z                                                              )                                                          (        3        )                                Equation        ⁢                                   ⁢        of        ⁢                                   ⁢        thermal        ⁢                                   ⁢        energy                                                                                   ∂                          θ              il                                            ∂            t                          =                                            -              u                        ⁢                                          ∂                                  θ                  il                                                            ∂                x                                              -                      v            ⁢                                          ∂                                  θ                  il                                                            ∂                y                                              +                                    ∂                              θ                il                                                    ∂              z                                +                                    ∂                              ∂                x                                      ⁢                          (                                                K                  h                                ⁢                                                      ∂                                          θ                      il                                                                            ∂                    x                                                              )                                +                                    ∂                              ∂                y                                      ⁢                          (                                                K                  h                                ⁢                                                      ∂                                          θ                      il                                                                            ∂                    y                                                              )                                +                                    ∂                              ∂                z                                      ⁢                          (                                                K                  h                                ⁢                                                      ∂                                          θ                      il                                                                            ∂                    t                                                              )                                +                                    (                                                ∂                                      θ                    il                                                                    ∂                  t                                            )                        rad                                              (        4        )                                Equation        ⁢                                   ⁢        of        ⁢                                   ⁢        moisture        ⁢                                   ⁢        diffusion                                                                                   ∂                          r              n                                            ∂            t                          =                                            -              u                        ⁢                                          ∂                                  r                  n                                                            ∂                x                                              -                      v            ⁢                                          ∂                                  r                  n                                                            ∂                y                                              -                      w            ⁢                                          ∂                                  r                  n                                                            ∂                z                                              +                                    ∂                              ∂                x                                      ⁢                          (                                                K                  h                                ⁢                                                      ∂                                          r                      n                                                                            ∂                    x                                                              )                                +                                    ∂                              ∂                y                                      ⁢                          (                                                K                  h                                ⁢                                                      ∂                                          r                      n                                                                            ∂                    y                                                              )                                +                                    ∂                              ∂                z                                      ⁢                          (                                                K                  h                                ⁢                                                      ∂                                          r                      n                                                                            ∂                    z                                                              )                                                          (        5        )                                Equation        ⁢                                   ⁢        of        ⁢                                   ⁢        continuity                                                                                        ⁢                                            ∂                              π                ′                                                    ∂              t                                =                                    -                                                R                  ⁢                                                                           ⁢                                      π                    0                                                                                        c                    v                                    ⁢                                      ρ                    0                                    ⁢                                      θ                    0                                                                        ⁢                          (                                                                                          ∂                                              ρ                        0                                                              ⁢                                          θ                      0                                        ⁢                    u                                                        ∂                    x                                                  +                                                                            ∂                                              ρ                        0                                                              ⁢                                          θ                      0                                        ⁢                    v                                                        ∂                    y                                                  +                                                                            ∂                                              ρ                        0                                                              ⁢                                          θ                      0                                        ⁢                    w                                                        ∂                    z                                                              )                                                          (        6        )            where    u, v, w: wind speeds    f: Corioli's parameter    Km: eddy viscosity coefficient of momentum    Kh: eddy diffusion coefficient of heat and moisture    θil: potential temperature of moisture (ice-water)    rn: mixing ratio of moisture such as rain or snow    ρ: density    rad: radiation    g: gravitational acceleration    π′: Exner function (change)    θv: temporary potential temperature    p: pressure    R: gas constant    Cv: isovolumic specific heat    Subscript 0: reference value
In this manner, computations are made according to the fundamental equations for wind field analysis represented by the RAMS (Regional Atmospheric Modeling System) code to obtain wind direction vector data (wind field data) showing the meteorological data at the respective parent grid point positions and the meteorological data at the respective child grid point positions at predetermined time intervals (for example, 20-second intervals) from the start of computation.
<Explanation for Outline of Diffusion Status Predictive Computation>
Next, the diffusion status predictive computation will be explained. To perform the diffusion status predictive computation, the wind field data at the respective parent grid point positions and the respective child grid point positions at 20-second intervals, which were obtained by the RAMS (Regional Atmospheric Modeling System) code, were substituted, one after another, into the HYPACT (Hybrid Particle Concentration Transport Model) code developed by Colorado State University and ATMET. As a concrete example of the diffusion status predictive computation, Lagrangian particle dispersion model is adopted.
In the Lagrangian particle dispersion model, the diffusion velocities (u′, v′ and w′) are calculated using equations (7) to (9) indicated below, and the respective particles are moved.u′(t)=Ruu′(t−Δt)+(1−Ru2)ruv′(t)=Rvv′(t−Δt)+(1−Rv2)rvw′(t)=Rww′(t−Δt)+(1−Rw2)rw  (7)where Ru, Rv and Rw: Lagrange turbulent autocorrelation coefficient,
u′(t), v′(t) and w′(t): turbulent diffusion velocity component of particle, and
t: time.                                                         R              u                        ⁡                          (                              Δ                ⁢                                                                   ⁢                t                            )                                =                                                                                                                u                      ′                                        ⁡                                          (                      t                      )                                                        ·                                                            u                      ′                                        ⁡                                          (                                              t                        -                                                  Δ                          ⁢                                                                                                           ⁢                          t                                                                    )                                                                      _                                            σ                u                2                                      =                          exp              ⁡                              (                                  -                                                            Δ                      ⁢                                                                                           ⁢                      t                                                              T                      Lu                                                                      )                                                    ⁢                                  ⁢                                            R              v                        ⁡                          (                              Δ                ⁢                                                                   ⁢                t                            )                                =                                                                                                                v                      ′                                        ⁡                                          (                      t                      )                                                        ·                                                            v                      ′                                        ⁡                                          (                                              t                        -                                                  Δ                          ⁢                                                                                                           ⁢                          t                                                                    )                                                                      _                                            σ                v                2                                      =                          exp              ⁡                              (                                  -                                                            Δ                      ⁢                                                                                           ⁢                      t                                                              T                      Lv                                                                      )                                                    ⁢                                  ⁢                                            R              w                        ⁡                          (                              Δ                ⁢                                                                   ⁢                t                            )                                =                                                                                                                w                      ′                                        ⁡                                          (                      t                      )                                                        ·                                                            w                      ′                                        ⁡                                          (                                              t                        -                                                  Δ                          ⁢                                                                                                           ⁢                          t                                                                    )                                                                      _                                            σ                w                2                                      =                          exp              ⁡                              (                                  -                                                            Δ                      ⁢                                                                                           ⁢                      t                                                              T                      Lw                                                                      )                                                                        (        8        )            where σu, σv, σw: turbulent velocity standard deviation,
TLu, TLv, TLw: Lagrange time scale,ru=σuηu, rv=σvηv, rw=σwηw+wd  (9)where ηu, ηv, ηw: normal random number (the mean is zero),
wd: gravitational sedimentation rate.
A concrete example will be described in which the wind field data at the respective parent grid point positions and the respective child arid point positions at 20-second intervals, which had been obtained by the RAMS (Regional Atmospheric Modeling System) code, were substituted, one after another, into the HYPACT (Hybrid Particle Concentration Transport Model) code to perform the diffusion status predictive computation.
To carry out this computation, a setting is made such that the substance emitted from the emission source into the atmosphere is converted into many particles P, and N (for example, twenty) particles P are generated from the position of the emission source in each computation cycle of Δt (here Δt=20 seconds).
Concretely, 20 particles P are generated in each computation cycle of Δt (20 seconds) such that 20 particles are generated at the start of computation, 20 particles are generated 20 seconds after the start of computation, 20 particles are generated 40 seconds after the start of computation, and so on. Also, the positions (spatial coordinates) of the respective particles P are found by computation in each computation cycle of Δt (20 seconds).
The 20 particles P generated at the start of computation (time 0) are designated as P0001, P0002, P0003, P0004, P0005, P0006, P0007, P0008, P0009, P0010, P0011, P0012, P0013, P0014, P0015, P0016, P0017, P0018, P0019, and P0020.
The 20 particles P generated 20 seconds after the start of computation are designated as P2001, P2002, P2003, P2004, P2005, P2006, P2007, P2008, P2009, P2010, P2011, P2012, P2013, P2014, P2015, P2016, P2017, P2018, P2019, and P2020.
The 20 particles P generated 40 seconds after the start of computation are designated as P4001, P4002, P4003, P4004, P4005, P4006, P4007, P4008, P4009, P4010, P4011, P4012, P4013, P4014, P4015, P4016, P4017, P4018, P4019, P4020.
In the above designations, the figure shown in the lower row behind the symbol “P” represents the point of time elapsed after the start of computation, while the figure shown in the upper row behind the symbol “P” represents the particular particle of the twenty particles generated at this point of time. The particles generated at other points of time are designated similarly.
In detail, at the start of computation, the following 20 particles are generated from the emission source S: P0001, P0002, P0003, P0004, P0005, P0006, P0007, P0008, P0009, P0010, P0011, P0012, P0013, P0014, P0015, P0016, P0017, P0018, P0019, and P0020.
Twenty seconds after the start of computation, the following 20 particles are newly generated from the emission source S shown in FIG. 19: P2001, P2002, P2003, P2004, P2005, P2006, P2007, P2008, P2009, P2010, P2011, P2012, P2013, P2014, P2015, P2016, P2017, P2018, P2019, and P2020.
At this time, the particles generated at the start of computation, i.e., P0001, P0002, P0003, P0004, P0005, P0006, P0007, P0008, P0009, P0010, P0011, P0012, P0013, P0014, P0015, P0016, P0017, P0018, P0019, and P0020, arrive at positions remote from the emission source S while diffusing.
The positions of the respective particles P are determined by calculating the diffusion velocities (u′, v′, w′) of the respective particles P in the Lagrangian particle dispersion model with the use of the wind field data at 20-second intervals obtained by the RAMS (Regional AtmosphericModeling System) code, and moving the respective particles based on the results.
Forty seconds after the start of computation, the following 20 particles are newly generated from the emission source S shown in FIG. 20: P4001, P4002, P4003, P4004, P4005, P4006, P4007, P4008, P4009, P4010, P4011, P4012, P4013, P4014, P4015, P4016, P4017, P4018, P4019, P4020.
At this time, the particles generated at the start of computation, i.e., P0001, P0002, P0003, P0004, P0005, P0006, P0007, P0008, P0009, P0010, P0011, P0012, P0013, P0014, P0015, P0016, P0017, P0018, P0019, and P0020, reach positions further remote from the emission source S while diffusing further.
Moreover, the following 20 particles generated 20 seconds after the start of computation, P2001, P2002, P2003, P2004, P2005, P2006, P2007, P2008, P2009, P2010, P2011, P2012, P2013, P2014, P2015, P2016, P2017, P2018, P2019, and P2020, arrive at positions remote from the emission source S while diffusing.
The positions of the respective particles P are determined by calculating the diffusion velocities (u′, v′, w′) of the respective particles P in the Lagrangian particle dispersion model with the use of the wind field data at 20-second intervals obtained by the RAMS (Regional Atmospheric Modeling System) code, and moving respective particles based on the results.
Sixty seconds after the start of computation, following 20 particles are newly generated from the emission source S shown in FIG. 21: P6001, P6002, P6003, P6004, P6005, P6006, P6007, P6008, P6009, P6010, P6011, P6012, P6013, P6014, P6015, P6016, P6017, P6018, P6019, and P6020.
At this time, the particles generated at the of computation, i.e. P0001, P0002, P0003, P0004, P0005, P0006, P0007, P0008, P0009, P0010, P0011, P0012, P0013, P0014, P0015, P0016, P0017, P0018, P0019, and P0020, reach positions further remote from the emission source S while diffusing further.
Moreover, the following 20 particles generated 20 seconds after the start of computation, P2001, P2002, P2003, P2004, P2005, P2006, P2007, P2008, P2009, P2010, P2011, P2012, P2013, P2014, P2015, P2016, P2017, P2018, P2019, and P2020, arrive at positions further remote from the emission source S while diffusing further.
Furthermore, the following 20 particles generated 40 seconds after the start of computation, P4001, P4002, P4003, P4004, P4005, P4006, P4007, P4008, P4009, P4010, P4011, P4012, P4013, P4014, P4015, P4016, P4017, P4018, P4019, P4020, reach positions remote from the emission source S while diffusing.
The positions of the respective particles P are determined by calculating the diffusion velocities (u′, v′, w′) of the respective particles P in the Lagrangian particle dispersion model with the use of the wind field data at 20-second intervals obtained by the RAMS (Regional Atmospheric Modeling System) code, and moving the respective particles based on the results.
As described above, groups of 20 particles are generated successively in respective computation cycles of Δt (20 seconds) and, at the same time, the position, i.e., spatial coordinates (xi(t), yi(t), zi(t)), of each particle in each computation cycle of Δt (20 seconds) is determined.
If the particles P are present in a unit space (a unit volume in a predicted area) separated by a predetermined distance from the emission source S when a predetermined period of time has elapsed since the start of computation, as shown in FIG. 22, the concentration of the substance in this unit space can be calculated from the number of these particles.
In detail, assume that the substance in an amount of Q (m3/second) is emitted at the emission source S. Since the particles P are generated at a rate of 20 particles in 20 seconds (namely, one particle per second), each particle P has an emission source intensity of Q/1 (m3). Thus, the number of the particles P present in the unit space is multiplied by the emission source intensity of Q/1 (m3), whereby the concentration of the substance in this unit space can be determined.
The foregoing concrete example will be described generally as follows: A substance, such as a gaseous body, emitted from the emission source is replaced by many particles. N particles per second are emitted from the emission source. In this case, the calculated amount of emission of the particles is N/sec. If the actual amount of emission of the substance emitted from the emission source is Q (m3/sec), each particle has an emission source intensity of Q/N (m3).
The equations of particle motion are numerically calculated in an unsteady manner for each particle. That is, the wind field data obtained by the RAMS (Regional Atmosphere Modeling System) code are substituted into the HYPACT (Hybrid Particle Concentration Transport Model) code which covers the equations of motion of particles, the diffusion velocities (u′, v′, w′) of each particle P are calculated using the Lagrangian particle dispersion model, and the respective particles are moved. In this manner, the coordinates of each particle can be determined unsteadily. In other words, the spatial coordinates of each particle can be determined in each computation cycle of Δt. The data on each particle, obtained by the Lagrangian particle model and recorded into a data recorder, are only the spatial coordinates (xi(t), yi(t), zi(t)) of each particle.
The HYPACT code, which represents the equations of motion of particles (substance), expresses the advection, diffusion and gravitational sedimentation phenomena of particles. The advection phenomenon of particles depends on the time-averaged velocity of the atmosphere, the diffusion phenomenon of particles depends on the turbulent velocity of the atmosphere, and the gravitational sedimentation of particles depends on the mass of particles, gravitational acceleration, and viscosity coefficient of air (see FIG. 23).
If the number of particles in the unit volume in the air is n, the gas concentration (substance concentration) in this space is n×Q/N (gas M3/air m3), namely, the product of the number n of particles present in this unit space and the emission source intensity Q/N which each particle has.
This environmental concentration (the substance concentration in the unit volume) depends on timewise variations in the amount of emission of the substance emitted. Thus, under conditions where the amount of emission varies with time, calculation of diffusion needs to be made for each of the emission conditions. If there are many emission conditions which can be anticipated, it is necessary to carry out calculations for diffusion corresponding to different cases of the amount of emission. This requires a huge calculation time.
That is, as shown in FIG. 24, when a gas (substance) is emitted from an emission source S (e.g., a chimney), for example, timewise variations in the gas concentration at a location F downwind of the emission source S depend on timewise variations in the substance emitted from the emission source S.
In detail, if the amount of emission of the substance varies with time as in FIG. 25(a), the concentration of the substance at the location F varies with time as in FIG. 25(b). If the amount of emission of the substance is constant as in FIG. 26(a), the concentration of the substance at the location F rises to a constant value and then remains constant as in FIG. 26(b). If the substance is emitted instantaneously as in FIG. 27(a), the concentration of the substance at the location F rises transiently and then falls to zero as in FIG. 27(b).
If the amount of emission of the substance varies with the passage of time, as described above, the number of the particles generated needs to be changed over time in correspondence with the amount of emission of the substance. Then, the positions which the generated particles, whose number was changed with the passage of time in this manner, are moved to are found. From these positions to which the particles moved, the concentration of the substance is calculated. Thus, calculations of diffusion must be carried out for different cases where the amount of emission varies by itself and over time. This necessitates enormous calculation results.
In an accident in which a radioactive substance is emitted to the outside at facilities handling the radioactive substance, for example, a very large number of substances (e.g., about 100 kinds of substances) are emitted. Moreover, the amount of emission of each substance differs according to time. Thus, the number of particles generated is changed for each substance with the passage of time in correspondence with the amount of emission of each substance. A determination is made of where the generated particles, thus changed in number, are moved to. From these positions to which the particles are moved, the concentration of the substance is calculated. In case of the accident, therefore, there is need to carry out 100 kinds of diffusion calculations corresponding to 100 kinds of substances, for example.